Préparation

Packages

library(dplyr)
library(ggplot2)
library(plotly)
library(glmnet)
library(gbm)

Chargement des données

load("../1_data_management_dplyr/fichiers_prepared.RData")

Jointure INSEE & FINESS

nb_ET_CODGEO=finess_et%>%
  group_by(CODGEO)%>%
  summarise(nb_ET=n())

sum(!nb_ET_CODGEO$CODGEO%in%insee$CODGEO)
## [1] 3191
nb_ET_CODGEO=merge(nb_ET_CODGEO,insee,by="CODGEO",all.y=T)
nb_ET_CODGEO <- nb_ET_CODGEO%>%mutate(nb_ET=ifelse(is.na(nb_ET),0,nb_ET))
table(nb_ET_CODGEO$nb_ET)%>%head
## 
##     0     1     2     3     4     5 
## 28153  2888  1400   947   583   425

Sélection des variables pour entraîner un modèle

data_model <-  nb_ET_CODGEO%>%
  select(-LIBGEO,-CODGEO)%>%
  mutate_if(is.character,as.numeric)

Echantillonage pour séparer apprentissage/test ie train / test

train=sample(x = 1:nrow(data_model),size = round(.65*nrow(data_model)))

GLM

Préparation des données

Les valeurs manquantes sont très mal gérées par les GLM, on va donc imputer par la moyenne.
Il existe de nombreuses stratégies d’imputation donc certaines s’appuient sur des modèles de machine learning avancé : - Moyenne - Médiane - HotDeck - HotDeck stratifié - k plus proches voisins (kNN) - Régression/classification GLM - Régression/classification CART - Régression/classification GBM/RandomForest

Il existe plusieurs packages R qui implémente des stratégies d’imputation comme MICE ou simputation.
Lorsque les contraintes et le contexte sont très spécifiques et demande une maîtrise importante, il ne faut pas hésiter à implémenter son modèle d’imputation.
Vous pouvez vous référer à la présentation (Séminaire EIG #2) sur l’imputation des non-réponses partielles pour l’enquête OC : https://gitlab.com/DREES_code/OSAM/enquete_oc

data_model_imputed_for_glm=data_model%>%
  mutate_all(function(x)ifelse(is.na(x),mean(x,na.rm=T),x))

Entraînement du modèle

Un GLM dans R ce décrit comme une formule Y~X1+X2+X3 où Y est la variable à expliquer/prédire et X1, X2, X4 les variables explicatives.
Lorsqu’on a déjà sélectionné les variables pertinentes pour le modèle, il suffit d’écrit . pour utiliser toutes les variables disponibles sauf Y.
Si on a souhaité conserver la variable ID et qu’on souhaite garder toutes les variables sauf ID il suffit d’écrire Y~.-ID.

On peut ajuster le modèle avec un certain nombre de parmaètres facultatifs. - Choisir la fonction de lien (log, inverse, logit, identité…) et la loi conditionnelle. Si on a oublié quelle est la fonction de lien canonique associée à une loi de la famille exponentielle, pas de problème, elle est assignée par défaut. - On peut préciser un ou plusieurs offset (coeff fixé, pas forcément à 1). - Pondérer les observations avec weights (par exemple si on veut donner moins d’importance aux données plus anciennes). Des idées générales sur cette approche, valable pour la plupart des modèles prédictifs. ici. - On peut définir les contrasts pour préciser quelle catégorie doit être utilisée comme référence pour une variable explicative catégorielle.

model <- glm(nb_ET~.,
             data=data_model_imputed_for_glm[train,],
             family=poisson(link = "log"))

Coefficients du modèle

coeff_glm=summary(model)$coefficients
coeff_glm
##                      Estimate   Std. Error    z value      Pr(>|z|)
## (Intercept)     -1.126276e+01 1.052242e+00 -10.703581  9.791915e-27
## NBMENFISC14     -5.527830e-05 1.609641e-06 -34.342000 1.854201e-258
## NBPERSMENFISC14  2.694129e-05 7.808639e-07  34.501902 7.511508e-261
## MED14            3.515499e-05 2.674837e-06  13.142852  1.870430e-39
## PIMP14           1.081355e-01 3.319893e-03  32.571972 1.023125e-232
## TP6014          -3.693168e-01 6.034702e-03 -61.198850  0.000000e+00
## TP60AGE114      -5.709656e-03 1.773697e-03  -3.219070  1.286070e-03
## TP60AGE214      -5.125060e-02 2.493470e-03 -20.553926  7.098366e-94
## TP60AGE314      -6.478064e-03 2.846824e-03  -2.275541  2.287349e-02
## TP60AGE414       1.785775e-02 3.142301e-03   5.683017  1.323390e-08
## TP60AGE514       2.139041e-02 3.350546e-03   6.384157  1.723447e-10
## TP60AGE614      -1.435137e-01 2.827824e-03 -50.750585  0.000000e+00
## TP60TOL114      -1.857351e-02 4.571160e-03  -4.063195  4.840554e-05
## TP60TOL214      -5.768398e-02 2.596958e-03 -22.212134 2.621841e-109
## PRA14            1.862361e-01 7.775149e-03  23.952742 8.651937e-127
## PTSAC14         -1.303710e-02 5.698780e-03  -2.287701  2.215494e-02
## PPEN14           2.039243e-01 9.969454e-03  20.454914  5.431799e-93
## PPAT14           2.629440e-01 9.632152e-03  27.298574 4.410173e-164
## PPSOC14         -2.113385e-01 1.078465e-01  -1.959624  5.003979e-02
## PPFAM14          7.672117e-01 1.076224e-01   7.128736  1.012950e-12
## PPMINI14         5.105354e-01 1.095562e-01   4.660033  3.161587e-06
## PPLOGT14         3.224468e+00 1.121298e-01  28.756559 7.499055e-182
## D114            -1.010405e-03 2.484651e-05 -40.665886  0.000000e+00
## D914            -1.485463e-04 5.893037e-06 -25.207092 3.349085e-140
## RD14             1.362044e+00 6.057189e-02  22.486407 5.638491e-112

D’autres visualisations 

# install.packages("stargazer")
stargazer::stargazer(model, type = "html")
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>nb_ET</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">NBMENFISC14</td><td>-0.0001<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">NBPERSMENFISC14</td><td>0.00003<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">MED14</td><td>0.00004<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PIMP14</td><td>0.108<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP6014</td><td>-0.369<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.006)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE114</td><td>-0.006<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.002)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE214</td><td>-0.051<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.002)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE314</td><td>-0.006<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE414</td><td>0.018<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE514</td><td>0.021<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE614</td><td>-0.144<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60TOL114</td><td>-0.019<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.005)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60TOL214</td><td>-0.058<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PRA14</td><td>0.186<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.008)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PTSAC14</td><td>-0.013<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.006)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PBEN14</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPEN14</td><td>0.204<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.010)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPAT14</td><td>0.263<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.010)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPSOC14</td><td>-0.211<sup>*</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.108)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPFAM14</td><td>0.767<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.108)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPMINI14</td><td>0.511<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.110)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPLOGT14</td><td>3.224<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.112)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PIMPOT14</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">D114</td><td>-0.001<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00002)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">D914</td><td>-0.0001<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">RD14</td><td>1.362<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.061)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-11.263<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.052)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>23,804</td></tr>
## <tr><td style="text-align:left">Log Likelihood</td><td>-60,696.850</td></tr>
## <tr><td style="text-align:left">Akaike Inf. Crit.</td><td>121,443.700</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Prediction

pred_glmtrain=data.frame(pred=predict(model,newdata = data_model_imputed_for_glm[train,],type="response"), obs=data_model[train,]$nb_ET)
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading
perf_train=data.frame(pred=data_model[train,]$nb_ET, obs=data_model[train,]$nb_ET)


pred_glmtest=data.frame(pred=predict(model,newdata = data_model_imputed_for_glm[-train,],type="response"),obs=data_model[-train,]$nb_ET)
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading
perf_test=data.frame(pred=data_model[-train,]$nb_ET, obs=data_model[-train,]$nb_ET)

Courbe de Lorenz ou Gain Curves

Il s’agit des courbes de mesure des inégalités de richesse.
Ici on les applique à la mesure des inégalités de Y (dans notre cas nb_ET) dans la population : - La courbe “perfection” décrit les vraies d’inégalités - La courbe “random” décrit un modèle incapable de discerner les inégalités - La courbe “glm” décrit la capacité de notre modèle à appréhender (car échantillon de test) les inégalités

L’intérêt d’une telle mesure est qu’on peut définir - ce qu’est un “mauvais” modèle (random) - ce qu’est un “bon” modèle (perfection) et comparer notre modèle à ces deux extrêmes.
Cette approche est inspirée de l’étude de la courbe de ROC pour les modèles binaires (Bernoulli).

L’inconvénient est que cette métrique mesure les inégalités en rang et pas en taille.

pred_glmtest%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Lorenz_glm_points
  
pred_glmtest%>%
    na.omit%>%
    arrange(-obs)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Perfection
  
g <- ggplot()+
    geom_line(data = Lorenz_glm_points,
              aes(x=x100,y=y100,color="glm"))+
    geom_line(data = Perfection,
              aes(x=x100,y=y100,color="perfection")) +
    geom_line(data=data.frame(x100=c(0,1), yrandom=c(0,1)),
              aes(x=x100,y=yrandom,color="random"))
  
g

Coeff de Gini

Par extension de l’AUC sous la ROC, on calcule l’AUC sous la courbe de Lorenz.
Cette métrique, translatée (X->2*X-1) pour se comparer à “random” est appelée indice de Gini.
Pour bien faire, on peut rapporter le Gini du modèle au Gini de la “perfection” pour que notre indice ait des valeurs POSSIBLES entre 0 et 1.

Gini=function(pred){
  pred%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      # print(paste0("nombre_de_lignes:",nrow(.)," indice de Gini:"))
      2*mean(.$y)-1
      }
}
Gini(pred_glmtest)
## [1] 0.5972249
Gini(perf_test)
## [1] 0.9289132
Gini(pred_glmtrain)
## [1] 0.6228157
Gini(perf_train)
## [1] 0.9325213

Performance faible mais pas d’overfitting ie peu d’écart entre apprentissage et test.

GLMNET ou Elastic Net (mix L1 & L2)

Qu’est-ce que c’est ?

Très bonne explication ici. 2 paramètres à choisir : - alpha qui équilibre entre les pénalisations L1 et L2. - lambda qui équilibre la fonction de perte entre vraisemblance et pénalisation sur la taille des coefficients Lambda élevé signifie une pénalisation importante sur la taille des coefficients.
Lambda faible signifie qu’on est proche d’un GLM classique non pénalisé.

Le package R glmnet test automatiquement et de manière optimisé plusieurs lambda.
En revanche c’est à nous de tester plusieurs paramètres alpha. C’est ce qu’on appelle l’optimisation des hyperparamètres (hyperparameters tuning). La méthode qui consiste à tester les combinaisons des différents hyperparamètres s’appelle gridsearch.
La méthode gourmande/gloutonne (bruteforce) consiste à tester toutes les combinaisons.
Il existe des méthodes bayesiennes pour explorer intelligemment la grille de valeurs.
De nombreux packages R sont disponibles, mlrMBO est plutôt bien maintenu.
Mais ces considérations dépassent l’ambition de cette formation.

Préparation des données

Dans R il y a beaucoup de liberté quant au traitement des données avec les formats matrix, data.frame, data.table…
Pour glm, on avait un modèle de la forme cible ~ variables explicatives ou Y~X avec Y et X dans un même data.frame.
Pour glmnet, on doit fournir Y comme un vecteur et X comme une matrice de valeurs numériques.
Par conséquent les variables catégorielles devront être binarisées (dummy variables), ce procédé s’appelle le one-hot encoding. C’est très simple en R, il suffit d’utiliser la fonction model.matrix. Un exemple ici
Pas besoin d’y recourir dans notre cas. Si

target_train = data_model_imputed_for_glm[train,]$nb_ET
explanatory_train = data_model_imputed_for_glm[train,]%>%
  select(-nb_ET)%>%as.matrix

target_test = data_model_imputed_for_glm[-train,]$nb_ET
explanatory_test = data_model_imputed_for_glm[-train,]%>%
  select(-nb_ET)%>%as.matrix

Entraînement du modèle

thresh est le coefficient de convergence. On peut lire dans la documentation :
thresh Convergence threshold for coordinate descent. Each inner coordinate-descent
loop continues until the maximum change in the objective after any coefficient
update is less than thresh times the null deviance. Defaults value is 1E-7.

glmnet_model <- glmnet(x=explanatory_train,
       y=target_train,
       family="poisson",
       alpha=1,
       nlambda=100,#par défaut
       thresh = 1e-06,#par défaut
       maxit=1e7)
plot(glmnet_model$lambda,glmnet_model$dev.ratio)

# Distribution des coefficients en fonction du lambda.
plot(glmnet_model, xvar='lambda',label=T)

s=sample(glmnet_model$lambda,1)
smin=min(glmnet_model$lambda)
  pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=0)[,1],obs=target_test)
  print(paste("s=0",Gini(pred_glmnettest)))
## [1] "s=0 0.591052172672043"
  pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
  print(paste("s=smin=",smin,"Gini=",Gini(pred_glmnettest)))
## [1] "s=smin= 0.00156900416023852 Gini= 0.591052172672043"
    pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=s)[,1],obs=target_test)
  print(paste("s=random",s,"Gini=",Gini(pred_glmnettest)))
## [1] "s=random 0.00249829832594608 Gini= 0.590321105877747"

[Avancé] cross-validation et choix du lambda optimal

Pour automatiser et rationnaliser le choix du lambda on fait de la validation croisée (cross-validation) ie on coupe le dataset en NFOLDS, disons en 10, et entraîne sur 90% vs validation sur 10% 10 fois.

NFOLDS=10
for (alpha in 0:10/10){
  print(alpha)
  
  tryCatch(rm(glmnet_cv,pred_glmnettest),warning = function(w) {
    print(paste("warning",w))
}, error = function(e) {
    print(paste("error",e))
})
  tryCatch({
glmnet_cv = cv.glmnet(x = explanatory_train, y = target_train, 
                              family = "poisson",#'gaussian', 
                              # Pénalisation L1 100%
                              alpha = alpha,
                              # On s'intéresse à la deviance - on suppose que la distribution conditionnelle suit une loi de Poisson
                              type.measure = "deviance",#'rmse'
                              # 5-fold cross-validation
                              nfolds = NFOLDS,
                              # valeurs élevée pour un entraînement plus rapide mais moins performant
                              thresh = 1e-3,
                              # On limite le nombre d'itérations
                              maxit = 1e5)
s1se=glmnet_cv$lambda.1se
smin=glmnet_cv$lambda.min
  pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=0)[,1],obs=target_test)
  print(paste("s=0",Gini(pred_glmnettest)))
  pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
  print(paste("s=smin",Gini(pred_glmnettest)))
    pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=s1se)[,1],obs=target_test)
  print(paste("s=s1se",Gini(pred_glmnettest)))
  },warning = function(w) {
    print(paste("warning",w))
}, error = function(e) {
    print(paste("error",e))
})
}
## [1] 0
## [1] "warning simpleWarning in rm(glmnet_cv, pred_glmnettest): objet 'glmnet_cv' introuvable\n"
## [1] "s=0 0.530549879254173"
## [1] "s=smin 0.491033845212646"
## [1] "s=s1se -0.125987883560982"
## [1] 0.1
## [1] "s=0 0.569958405168113"
## [1] "s=smin 0.567633591437849"
## [1] "s=s1se 0.543620369155664"
## [1] 0.2
## [1] "s=0 0.570660871328051"
## [1] "s=smin 0.56847105691447"
## [1] "s=s1se 0.556439987548278"
## [1] 0.3
## [1] "s=0 0.571030222945907"
## [1] "s=smin 0.569235536742053"
## [1] "s=s1se 0.547383521034902"
## [1] 0.4
## [1] "s=0 0.571227063029675"
## [1] "s=smin 0.569423306750592"
## [1] "s=s1se 0.546153643328189"
## [1] 0.5
## [1] "s=0 0.57165069821929"
## [1] "s=smin 0.569992293677941"
## [1] "s=s1se 0.543537715483627"
## [1] 0.6
## [1] "s=0 0.571429104538633"
## [1] "s=smin 0.569558630255826"
## [1] "s=s1se 0.594412732800699"
## [1] 0.7
## [1] "s=0 0.571385886602735"
## [1] "s=smin 0.569975162114087"
## [1] "s=s1se 0.546906481184709"
## [1] 0.8
## [1] "s=0 0.571488474268542"
## [1] "s=smin 0.569781563915826"
## [1] "s=s1se 0.551642407637421"
## [1] 0.9
## [1] "s=0 0.571581127356319"
## [1] "s=smin 0.56966485603759"
## [1] "s=s1se 0.550196447455412"
## [1] 1
## [1] "s=0 0.572030488708477"
## [1] "s=smin 0.570142709961255"
## [1] "s=s1se 0.548101173627391"
Gini(pred_glmtest)
## [1] 0.5972249

Comparaison pour différents lambda

plot(glmnet_cv)

Coefficients du modèle pour le “meilleur” lambda

Des idées pratiques en anglais ici
et théoriques en français ici

s1se=glmnet_cv$lambda.1se
smin=glmnet_cv$lambda.min
coeff_glmnet=coef(glmnet_model,s=s1se)
coeff_glmnet
## 27 x 1 sparse Matrix of class "dgCMatrix"
##                             1
## (Intercept)      1.353915e+01
## NBMENFISC14     -4.239680e-08
## NBPERSMENFISC14  2.691878e-07
## MED14            1.551383e-05
## PIMP14           7.410487e-02
## TP6014          -2.098730e-01
## TP60AGE114       .           
## TP60AGE214      -4.625761e-02
## TP60AGE314      -4.476880e-04
## TP60AGE414       4.092152e-03
## TP60AGE514       8.465594e-03
## TP60AGE614      -1.363710e-01
## TP60TOL114      -4.754384e-02
## TP60TOL214      -7.678858e-02
## PRA14           -1.093402e-02
## PTSAC14          .           
## PBEN14           .           
## PPEN14           .           
## PPAT14           2.277342e-02
## PPSOC14          .           
## PPFAM14          7.653354e-02
## PPMINI14         .           
## PPLOGT14         2.198343e+00
## PIMPOT14        -1.571213e-01
## D114            -1.243651e-03
## D914            -1.234753e-05
## RD14             .

Comparaison des modèles GLM vs GLMNET

Comparaison des coefficients

# On prépare la table des coeffs de glmnet
coeff_glmnet=as.matrix(coeff_glmnet)
coeff_glmnet=data.frame(var=rownames(coeff_glmnet),coeff_glmnet=coeff_glmnet[,1])
coeff_glmnet$var=as.character(coeff_glmnet$var)

# On prépare la table des coeffs de glm
coeff_glm=data.frame(var=rownames(coeff_glm),coeff_glm=coeff_glm[,1])
coeff_glm$var=as.character(coeff_glm$var)

# On joint les deux tables et on les compare
comparaison_coeff=merge(coeff_glm,coeff_glmnet,by="var")
comparaison_coeff%>%
  mutate(ratio_glm_sur_glmnet=coeff_glm/coeff_glmnet)%>%
  as.tbl
## # A tibble: 25 x 4
##    var               coeff_glm coeff_glmnet ratio_glm_sur_glmnet
##    <chr>                 <dbl>        <dbl>                <dbl>
##  1 (Intercept)     -11.3            1.35e+1               -0.832
##  2 D114             -0.00101       -1.24e-3                0.812
##  3 D914             -0.000149      -1.23e-5               12.0  
##  4 MED14             0.0000352      1.55e-5                2.27 
##  5 NBMENFISC14      -0.0000553     -4.24e-8             1304.   
##  6 NBPERSMENFISC14   0.0000269      2.69e-7              100.   
##  7 PIMP14            0.108          7.41e-2                1.46 
##  8 PPAT14            0.263          2.28e-2               11.5  
##  9 PPEN14            0.204          0.                   Inf    
## 10 PPFAM14           0.767          7.65e-2               10.0  
## # ... with 15 more rows

Corrélation entre les variables explicatives

Pour bien comprendre ce qui se passe, on jette un oeil à la matrice des corrélations. On se rend compte que lorsque deux variables sont très corrélées, en général le glmnet n’apprend que sur l’une des deux ie le coeff de l’autre passe à 0 ! C’est de la sélection de variables

cor(explanatory_train)%>%as.data.frame->cormat
rownames(cormat) <- colnames(cormat)
knitr::kable(cormat)
NBMENFISC14 NBPERSMENFISC14 MED14 PIMP14 TP6014 TP60AGE114 TP60AGE214 TP60AGE314 TP60AGE414 TP60AGE514 TP60AGE614 TP60TOL114 TP60TOL214 PRA14 PTSAC14 PBEN14 PPEN14 PPAT14 PPSOC14 PPFAM14 PPMINI14 PPLOGT14 PIMPOT14 D114 D914 RD14
NBMENFISC14 1.0000000 0.9980597 0.0215953 0.0028921 0.1052056 -0.0860245 -0.0521292 0.0070649 0.0158168 -0.0142194 -0.1186299 -0.0522551 0.0150650 0.0525470 0.0455294 0.0217454 -0.0628268 0.0333159 0.0691816 -0.0199791 0.0866507 0.1057856 -0.1308763 -0.1385921 0.1068657 0.2617827
NBPERSMENFISC14 0.9980597 1.0000000 0.0229011 0.0019599 0.1172289 -0.0808335 -0.0463743 0.0150583 0.0273077 0.0018635 -0.1046563 -0.0389131 0.0225112 0.0648001 0.0592228 0.0118646 -0.0786364 0.0229949 0.0823872 -0.0010762 0.0958551 0.1146606 -0.1300055 -0.1458092 0.1049247 0.2665857
MED14 0.0215953 0.0229011 1.0000000 0.4077471 -0.2910832 -0.1500026 -0.1992152 -0.2260350 -0.1888196 -0.1312729 -0.0769175 -0.1583395 -0.2742843 0.2530778 0.2292955 0.0560573 -0.1551675 0.1621254 -0.3739791 -0.3134540 -0.3341195 -0.3606243 -0.3806792 0.3803337 0.4262154 0.1411290
PIMP14 0.0028921 0.0019599 0.4077471 1.0000000 -0.6622976 -0.3469204 -0.4525352 -0.5028118 -0.4179267 -0.2907931 -0.1771124 -0.3892491 -0.6180156 0.5867046 0.5423128 0.0777533 -0.2849247 0.1704712 -0.8006538 -0.5981922 -0.7543827 -0.7850175 -0.8155333 0.8272099 0.6691432 0.0472803
TP6014 0.1052056 0.1172289 -0.2910832 -0.6622976 1.0000000 0.4284438 0.5973560 0.7048388 0.5996171 0.4650306 0.2822975 0.5854329 0.8107210 -0.3231671 -0.2910016 -0.0803133 0.0956991 -0.2073212 0.8369567 0.5302835 0.8542956 0.8167357 0.4226966 -0.7394627 -0.3623792 0.2468963
TP60AGE114 -0.0860245 -0.0808335 -0.1500026 -0.3469204 0.4284438 1.0000000 0.6996132 0.5802272 0.6210235 0.4895079 0.3840939 0.4727018 0.4755313 -0.2137757 -0.2080100 0.0222543 0.0734345 -0.1167832 0.4565836 0.3471938 0.4583764 0.3999994 0.3083101 -0.2192054 -0.2596827 -0.0823605
TP60AGE214 -0.0521292 -0.0463743 -0.1992152 -0.4525352 0.5973560 0.6996132 1.0000000 0.8224020 0.7893669 0.5453349 0.3082232 0.5997015 0.6355497 -0.2753687 -0.2578959 -0.0201535 0.1056694 -0.1578076 0.5852137 0.4322669 0.5901522 0.5218786 0.3844230 -0.3368280 -0.3315089 -0.0611591
TP60AGE314 0.0070649 0.0150583 -0.2260350 -0.5028118 0.7048388 0.5802272 0.8224020 1.0000000 0.7907886 0.5511260 0.3137646 0.6514217 0.7346698 -0.2727765 -0.2521532 -0.0360734 0.0860176 -0.1825043 0.6461682 0.4656216 0.6521336 0.5870576 0.4062684 -0.4339575 -0.3475086 0.0125118
TP60AGE414 0.0158168 0.0273077 -0.1888196 -0.4179267 0.5996171 0.6210235 0.7893669 0.7907886 1.0000000 0.7307591 0.4423994 0.7224931 0.5913820 -0.1538919 -0.1390782 -0.0357988 -0.0263840 -0.1879764 0.5598322 0.4517782 0.5675372 0.4582539 0.3626659 -0.3388408 -0.2962708 0.0043506
TP60AGE514 -0.0142194 0.0018635 -0.1312729 -0.2907931 0.4650306 0.4895079 0.5453349 0.5511260 0.7307591 1.0000000 0.6496638 0.7462134 0.4171335 -0.0159521 -0.0118301 -0.0162803 -0.1432056 -0.1691843 0.4257486 0.3814942 0.4536359 0.2806997 0.2724857 -0.2270215 -0.1962145 0.0273695
TP60AGE614 -0.1186299 -0.1046563 -0.0769175 -0.1771124 0.2822975 0.3840939 0.3082232 0.3137646 0.4423994 0.6496638 1.0000000 0.6533659 0.2502127 0.0232300 0.0131133 0.0437007 -0.1285273 -0.1326114 0.2720856 0.2327039 0.3504872 0.1039938 0.1872943 -0.1129029 -0.1122044 0.0167602
TP60TOL114 -0.0522551 -0.0389131 -0.1583395 -0.3892491 0.5854329 0.4727018 0.5997015 0.6514217 0.7224931 0.7462134 0.6533659 1.0000000 0.5554042 -0.1217988 -0.1263670 0.0508431 -0.0506555 -0.1472988 0.4968401 0.3935647 0.5602870 0.3325648 0.3090969 -0.3250373 -0.2132905 0.0812710
TP60TOL214 0.0150650 0.0225112 -0.2742843 -0.6180156 0.8107210 0.4755313 0.6355497 0.7346698 0.5913820 0.4171335 0.2502127 0.5554042 1.0000000 -0.3389176 -0.3123168 -0.0495672 0.1307905 -0.2078403 0.7202636 0.5322553 0.7134434 0.6614165 0.4911544 -0.5715571 -0.4330499 0.0210353
PRA14 0.0525470 0.0648001 0.2530778 0.5867046 -0.3231671 -0.2137757 -0.2753687 -0.2727765 -0.1538919 -0.0159521 0.0232300 -0.1217988 -0.3389176 1.0000000 0.9789818 -0.1330281 -0.8648094 -0.3531594 -0.3285947 0.0254452 -0.4175110 -0.4277412 -0.4454008 0.4861982 0.3883575 0.0518605
PTSAC14 0.0455294 0.0592228 0.2292955 0.5423128 -0.2910016 -0.2080100 -0.2578959 -0.2521532 -0.1390782 -0.0118301 0.0131133 -0.1263670 -0.3123168 0.9789818 1.0000000 -0.3323673 -0.8615756 -0.4084203 -0.2677540 0.0874546 -0.3730935 -0.3650827 -0.3557900 0.4478249 0.3207720 0.0099472
PBEN14 0.0217454 0.0118646 0.0560573 0.0777533 -0.0803133 0.0222543 -0.0201535 -0.0360734 -0.0357988 -0.0162803 0.0437007 0.0508431 -0.0495672 -0.1330281 -0.3323673 1.0000000 0.1876612 0.3516002 -0.2183883 -0.3073275 -0.1176672 -0.2039064 -0.3307316 0.0721418 0.2371115 0.1914874
PPEN14 -0.0628268 -0.0786364 -0.1551675 -0.2849247 0.0956991 0.0734345 0.1056694 0.0860176 -0.0263840 -0.1432056 -0.1285273 -0.0506555 0.1307905 -0.8648094 -0.8615756 0.1876612 1.0000000 0.1584739 0.0784269 -0.2469206 0.1906209 0.2033623 0.1828383 -0.2706725 -0.2480069 -0.0825585
PPAT14 0.0333159 0.0229949 0.1621254 0.1704712 -0.2073212 -0.1167832 -0.1578076 -0.1825043 -0.1879764 -0.1691843 -0.1326114 -0.1472988 -0.2078403 -0.3531594 -0.4084203 0.3516002 0.1584739 1.0000000 -0.3916229 -0.5438205 -0.2652198 -0.2918097 -0.3397489 0.1757270 0.4534294 0.3265344
PPSOC14 0.0691816 0.0823872 -0.3739791 -0.8006538 0.8369567 0.4565836 0.5852137 0.6461682 0.5598322 0.4257486 0.2720856 0.4968401 0.7202636 -0.3285947 -0.2677540 -0.2183883 0.0784269 -0.3916229 1.0000000 0.7775157 0.9418766 0.9502929 0.6102688 -0.8001881 -0.5919658 0.0197915
PPFAM14 -0.0199791 -0.0010762 -0.3134540 -0.5981922 0.5302835 0.3471938 0.4322669 0.4656216 0.4517782 0.3814942 0.2327039 0.3935647 0.5322553 0.0254452 0.0874546 -0.3073275 -0.2469206 -0.5438205 0.7775157 1.0000000 0.5513702 0.6205474 0.6262319 -0.4804263 -0.6514635 -0.3087437
PPMINI14 0.0866507 0.0958551 -0.3341195 -0.7543827 0.8542956 0.4583764 0.5901522 0.6521336 0.5675372 0.4536359 0.3504872 0.5602870 0.7134434 -0.4175110 -0.3730935 -0.1176672 0.1906209 -0.2652198 0.9418766 0.5513702 1.0000000 0.9076897 0.5034243 -0.7954744 -0.4576716 0.1678671
PPLOGT14 0.1057856 0.1146606 -0.3606243 -0.7850175 0.8167357 0.3999994 0.5218786 0.5870576 0.4582539 0.2806997 0.1039938 0.3325648 0.6614165 -0.4277412 -0.3650827 -0.2039064 0.2033623 -0.2918097 0.9502929 0.6205474 0.9076897 1.0000000 0.5383898 -0.8379558 -0.5243251 0.1148539
PIMPOT14 -0.1308763 -0.1300055 -0.3806792 -0.8155333 0.4226966 0.3083101 0.3844230 0.4062684 0.3626659 0.2724857 0.1872943 0.3090969 0.4911544 -0.4454008 -0.3557900 -0.3307316 0.1828383 -0.3397489 0.6102688 0.6262319 0.5034243 0.5383898 1.0000000 -0.5297794 -0.8423480 -0.4658049
D114 -0.1385921 -0.1458092 0.3803337 0.8272099 -0.7394627 -0.2192054 -0.3368280 -0.4339575 -0.3388408 -0.2270215 -0.1129029 -0.3250373 -0.5715571 0.4861982 0.4478249 0.0721418 -0.2706725 0.1757270 -0.8001881 -0.4804263 -0.7954744 -0.8379558 -0.5297794 1.0000000 0.5020017 -0.2934796
D914 0.1068657 0.1049247 0.4262154 0.6691432 -0.3623792 -0.2596827 -0.3315089 -0.3475086 -0.2962708 -0.1962145 -0.1122044 -0.2132905 -0.4330499 0.3883575 0.3207720 0.2371115 -0.2480069 0.4534294 -0.5919658 -0.6514635 -0.4576716 -0.5243251 -0.8423480 0.5020017 1.0000000 0.6621550
RD14 0.2617827 0.2665857 0.1411290 0.0472803 0.2468963 -0.0823605 -0.0611591 0.0125118 0.0043506 0.0273695 0.0167602 0.0812710 0.0210353 0.0518605 0.0099472 0.1914874 -0.0825585 0.3265344 0.0197915 -0.3087437 0.1678671 0.1148539 -0.4658049 -0.2934796 0.6621550 1.0000000

Prediction du GLMNET

pred_glmnettrain=data.frame(pred=predict(glmnet_model,newx = explanatory_train,type="response")[,1],obs=target_train)

pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)

Courbe de Lorenz

pred_glmnettest%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Lorenz_glmnet_points


g <- g +
    geom_line(data = Lorenz_glmnet_points,
              aes(x=x100,y=y100,color="glmnet"))
g%>%ggplotly

Coeff de Gini

Gini(pred_glmtest)
## [1] 0.5972249
Gini(pred_glmnettest)
## [1] 0.5910522
Gini(perf_test)
## [1] 0.9289132
Gini(pred_glmtrain)
## [1] 0.6228157
Gini(pred_glmnettrain)
## [1] 0.0150469
Gini(perf_train)
## [1] 0.9325213

Arbre de décision (regression et classification)

Commençons par utiliser un arbre de regression pour expliquer simplement le nombre d’établissements dans la commune. Ceci nous permettra de détecter les variables les plus influentes.

library(rpart)
options_rpart=rpart.control(minbucket =100, cp = 1E-3, maxdepth = 5)
tree <- rpart(nb_ET~., data=data_model_imputed_for_glm[train,],control = options_rpart)
rpart.plot::rpart.plot(tree,roundint=FALSE)

On offset la variable NBMENFISC14 pour calculer un nombre d’établissement par milliers de ménages. C’est aussi l’occasion d’introduire une autre visualisation plus interactive.

tree <- rpart(1000*nb_ET/NBMENFISC14~.-NBMENFISC14-NBPERSMENFISC14, data=data_model_imputed_for_glm[train,],control = options_rpart)

visNetwork::visTree(tree, main = "Arbre de décision interactif",
        width = "80%",  height = "400px")

GBM (gradient boosting method)

Préparation des données

L’avantage des modèles basés sur des arbres de décision, c’est qu’un arbre de décision sait bien gérer les NA.
Par exemple pour un arbre binaire, on peut imaginer séparer les NA dans un 3ème split, ou bien séparer NA vs le reste. La contrainte imposée par les variables continues est bien moins forte que dans un GLM. Avec un arbre si on coupe X à un seuil th (choisi par le modèle), on peut déciser de regrouper les NA avec le groupe 1 X<th ou avec le groupe 2 X>th.
L’imputation n’est pas nécessaire ici.

data_model%>%head
##   nb_ET NBMENFISC14 NBPERSMENFISC14   MED14  PIMP14  TP6014 TP60AGE114
## 1     0         304           799.5 21576.7      NA      NA         NA
## 2     0         100           235.5 21672.9      NA      NA         NA
## 3     0        6087         13660.5 19756.1 56.8215 15.7534    19.4181
## 4     0         603          1661.5 23204.8      NA      NA         NA
## 5     0          48           102.0 22157.5      NA      NA         NA
## 6     0        1030          2635.0 21679.3 61.7476  8.8425         NA
##   TP60AGE214 TP60AGE314 TP60AGE414 TP60AGE514 TP60AGE614 TP60TOL114
## 1         NA         NA         NA         NA         NA         NA
## 2         NA         NA         NA         NA         NA         NA
## 3    19.5204    19.1982    14.7159         NA         NA    5.40116
## 4         NA         NA         NA         NA         NA         NA
## 5         NA         NA         NA         NA         NA         NA
## 6         NA         NA         NA         NA         NA         NA
##   TP60TOL214 PRA14 PTSAC14 PBEN14 PPEN14 PPAT14 PPSOC14 PPFAM14 PPMINI14
## 1         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 2         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 3     24.796  71.8    67.9    3.9   27.3   10.1     6.5     2.8      1.8
## 4         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 5         NA    NA      NA     NA     NA     NA      NA      NA       NA
## 6         NA  76.3    73.8    2.5   26.5    8.2     4.3     2.7      0.8
##   PPLOGT14 PIMPOT14    D114    D914    RD14
## 1       NA       NA      NA      NA      NA
## 2       NA       NA      NA      NA      NA
## 3      1.9    -15.7 10545.7 33376.5 3.16495
## 4       NA       NA      NA      NA      NA
## 5       NA       NA      NA      NA      NA
## 6      0.8    -15.3 12797.2 34818.9 2.72082
params=c(shrinkage=.01,
         nb_trees=1000,
         depth=10,
         nminobs=10,
         bag.frac=.2)

gbm_model <- gbm(nb_ET~.
                 ,data=data_model[train,]
                 ,verbose=T
                 ,train.fraction=.8
                 ,shrinkage = params[1]
                 ,n.trees = params[2]
                 ,interaction.depth = params[3]
                 ,n.minobsinnode = params[4]
                 ,bag.fraction = params[5]
                 )
## Distribution not specified, assuming gaussian ...
## Iter   TrainDeviance   ValidDeviance   StepSize   Improve
##      1      161.6112        124.9445     0.0100    1.1945
##      2      159.5060        122.9226     0.0100    1.8978
##      3      158.3081        121.8718     0.0100    1.2360
##      4      156.5109        120.2592     0.0100    1.7279
##      5      154.6768        118.8279     0.0100    1.6574
##      6      152.7859        117.3525     0.0100    1.5123
##      7      151.5436        116.1283     0.0100    1.2773
##      8      150.0042        114.8612     0.0100    1.4503
##      9      147.9612        113.2020     0.0100    1.6624
##     10      146.6657        112.0619     0.0100    1.3079
##     20      135.3937        103.1946     0.0100    1.2281
##     40      115.0836         86.8909     0.0100    0.8974
##     60      101.1656         76.1005     0.0100    0.4093
##     80       91.7919         70.0738     0.0100    0.4968
##    100       84.9559         66.0688     0.0100    0.2363
##    120       80.7105         63.8572     0.0100    0.0554
##    140       76.8298         62.0454     0.0100    0.1331
##    160       72.9588         60.7204     0.0100   -0.1268
##    180       70.0501         60.2225     0.0100    0.0747
##    200       68.0310         60.3504     0.0100    0.1119
##    220       66.0278         60.5001     0.0100    0.0799
##    240       64.1627         60.8294     0.0100   -0.0563
##    260       63.4552         60.9274     0.0100   -0.0511
##    280       62.1450         61.4469     0.0100   -0.0137
##    300       61.1630         61.8902     0.0100   -0.0064
##    320       60.1808         62.1462     0.0100   -0.0569
##    340       59.3843         62.2452     0.0100    0.0246
##    360       58.6128         62.4991     0.0100   -0.0155
##    380       57.9869         62.9084     0.0100   -0.0024
##    400       57.0480         63.4925     0.0100    0.0210
##    420       56.5894         63.7226     0.0100    0.0317
##    440       55.9468         63.9451     0.0100   -0.0220
##    460       55.5446         63.9862     0.0100   -0.0205
##    480       54.9475         64.5412     0.0100   -0.0437
##    500       54.6453         64.6200     0.0100   -0.0058
##    520       54.2332         65.5485     0.0100   -0.0352
##    540       53.5622         66.1040     0.0100   -0.0117
##    560       52.8856         66.4188     0.0100    0.0406
##    580       52.2394         66.3664     0.0100    0.0133
##    600       51.7129         66.8059     0.0100   -0.0268
##    620       51.3457         66.8829     0.0100   -0.0309
##    640       50.8133         67.4678     0.0100   -0.0568
##    660       50.4482         67.2470     0.0100   -0.0526
##    680       49.9401         67.5536     0.0100    0.0125
##    700       49.4829         67.7143     0.0100   -0.0252
##    720       49.0172         68.7258     0.0100   -0.0031
##    740       48.5996         69.1680     0.0100   -0.0201
##    760       48.2110         69.2517     0.0100   -0.0419
##    780       47.8247         69.2563     0.0100   -0.0649
##    800       47.4762         69.2609     0.0100   -0.1056
##    820       47.0627         69.5084     0.0100   -0.0395
##    840       46.6427         69.5816     0.0100   -0.0627
##    860       46.4269         69.1736     0.0100   -0.0253
##    880       46.0515         69.3142     0.0100   -0.0526
##    900       45.8251         69.4874     0.0100   -0.0566
##    920       45.1928         69.9097     0.0100   -0.0063
##    940       44.9353         70.4944     0.0100   -0.0054
##    960       44.7112         70.7408     0.0100   -0.0065
##    980       44.4438         70.9160     0.0100    0.0046
##   1000       44.1010         71.2983     0.0100   -0.0346
pred_gbmtrain=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[train,]),
                        obs=target_train)
## Using 190 trees...
pred_gbmtest=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[-train,]),
                        obs=target_test)
## Using 190 trees...
print(paste("gbm test:",Gini(pred_gbmtest)))
## [1] "gbm test: 0.846249444647002"
print(paste("gbm train:",Gini(pred_gbmtrain)))
## [1] "gbm train: 0.862142145815112"
print(paste("glm:",Gini(pred_glmnettest)))
## [1] "glm: 0.591052172672043"

Il y a plus d’overfitting (écart de 1 point entre train et test) dans le GBM que dans le GLM, mais il y a aussi un fit de bien meilleur qualité, on va comparer les courbes GLM et GBM sur le test.

Courbe de Lorenz

pred_gbmtest%>%
    na.omit%>%
    arrange(-pred)%>%
    mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
    {
      data.frame(y100=quantile(.$y,0:100/100),
                 x100=quantile(.$x,0:100/100))
    }->Lorenz_gbm_points


g <- g +
    geom_line(data = Lorenz_gbm_points,
              aes(x=x100,y=y100,color="gbm"))
g%>%ggplotly

Interprétation d’un GBM

Importance des variables

vars=summary(gbm_model)$var%>%as.character
summary(gbm_model)

##                             var    rel.inf
## NBMENFISC14         NBMENFISC14 34.9627563
## NBPERSMENFISC14 NBPERSMENFISC14 18.8350015
## TP60AGE614           TP60AGE614 16.1073117
## TP60AGE114           TP60AGE114  4.9345014
## TP60AGE514           TP60AGE514  2.8233126
## RD14                       RD14  2.6426279
## TP60TOL114           TP60TOL114  1.6481053
## PPLOGT14               PPLOGT14  1.5172798
## PIMPOT14               PIMPOT14  1.5128370
## TP60AGE414           TP60AGE414  1.4323254
## PBEN14                   PBEN14  1.2352326
## D114                       D114  1.2160091
## TP60TOL214           TP60TOL214  1.1362209
## PPAT14                   PPAT14  1.1027992
## TP60AGE214           TP60AGE214  1.1016187
## PPEN14                   PPEN14  1.0699677
## TP60AGE314           TP60AGE314  1.0027726
## D914                       D914  0.8113241
## PRA14                     PRA14  0.8075889
## PTSAC14                 PTSAC14  0.7694872
## PPFAM14                 PPFAM14  0.6743689
## PPMINI14               PPMINI14  0.6011457
## TP6014                   TP6014  0.6000001
## PIMP14                   PIMP14  0.5668264
## MED14                     MED14  0.4791458
## PPSOC14                 PPSOC14  0.4094331

Effet marginal moyen

Attention, contrairement à un GLM, ces courbes ne s’interprètent pas comme des effets toutes choses égales par ailleurs parce qu’elles s’appuient sur la distribution réelle des données.

plot(gbm_model,i.var=vars[1])

plot(gbm_model,i.var=vars[2])

plot(gbm_model,i.var=vars[3])

plot(gbm_model,i.var=vars[4])

A comparer avec les valeurs réellement prises par ces variables, les extrêmes sont estimés mais souvent les valeurs sont aberrantes et non significatives.
## Monotonie des variables

Si nécessaire on peut imposer la monotonie des variables, ceci réduira le surapprentissage et donnera des effets plus propres/plus interprétables.

dans gbm : var.monotone il convient de donner un vecteur nommé (avec les noms des variables explicatives) et des valeurs -1 si décroissante, 1 si croissante, 0 si absence de contrainte.

monotony=rep(0,ncol(data_model)-1)
names(monotony) <- setdiff(names(data_model),"nb_ET")
monotony["TP60AGE114"] <- 1
params=c(shrinkage=.01,
         nb_trees=1000,
         depth=10,
         nminobs=10,
         bag.frac=.2)

gbm_model <- gbm(nb_ET~.
                 ,var.monotone = monotony
                 ,data=data_model[train,]
                 ,verbose=T
                 ,train.fraction=.8
                 ,shrinkage = params[1]
                 ,n.trees = params[2]
                 ,interaction.depth = params[3]
                 ,n.minobsinnode = params[4]
                 ,bag.fraction = params[5]
                 )
## Distribution not specified, assuming gaussian ...
## Iter   TrainDeviance   ValidDeviance   StepSize   Improve
##      1      161.4961        125.0582     0.0100    1.2216
##      2      159.8875        123.5936     0.0100    1.5888
##      3      158.4206        122.4782     0.0100    1.4064
##      4      157.5669        121.8906     0.0100    0.7681
##      5      155.6164        120.4929     0.0100    1.5742
##      6      154.0022        119.2011     0.0100    1.5202
##      7      152.0423        117.4708     0.0100    1.8606
##      8      151.0418        116.6187     0.0100    1.0506
##      9      149.9695        115.5649     0.0100    1.1009
##     10      148.6748        114.5025     0.0100    1.2641
##     20      135.6161        103.5849     0.0100    1.1795
##     40      115.8605         87.1145     0.0100    0.3968
##     60      101.1442         76.1938     0.0100    0.5126
##     80       91.3902         69.6246     0.0100    0.4239
##    100       85.1180         66.8352     0.0100    0.2402
##    120       80.5801         64.0562     0.0100   -0.0309
##    140       76.0434         62.3596     0.0100    0.3002
##    160       71.9982         61.3791     0.0100   -0.0483
##    180       69.7115         61.4357     0.0100    0.0863
##    200       67.8822         61.7476     0.0100    0.0993
##    220       66.3625         61.9588     0.0100   -0.0540
##    240       64.6558         62.1391     0.0100    0.0443
##    260       63.4095         62.4444     0.0100   -0.0294
##    280       62.4027         62.7356     0.0100   -0.0323
##    300       61.4847         63.2967     0.0100   -0.0543
##    320       60.1014         64.0942     0.0100   -0.0077
##    340       59.6189         64.2372     0.0100    0.0065
##    360       59.1495         64.5275     0.0100   -0.0536
##    380       58.2259         65.3635     0.0100   -0.0273
##    400       57.5878         65.7638     0.0100    0.0018
##    420       56.9295         66.4420     0.0100   -0.0161
##    440       56.3612         66.8646     0.0100   -0.0596
##    460       55.7183         67.6218     0.0100   -0.0187
##    480       55.3136         67.7511     0.0100   -0.0235
##    500       54.9858         67.5345     0.0100   -0.0296
##    520       54.3594         67.7049     0.0100    0.0438
##    540       53.9246         67.7988     0.0100   -0.0837
##    560       53.5223         68.1030     0.0100    0.0430
##    580       53.2301         67.7800     0.0100   -0.0834
##    600       52.8462         67.7087     0.0100   -0.0289
##    620       52.3831         67.9810     0.0100   -0.0006
##    640       51.9556         68.2019     0.0100   -0.0276
##    660       51.4675         68.1175     0.0100   -0.0459
##    680       51.2878         67.7366     0.0100    0.0050
##    700       50.7660         68.0417     0.0100   -0.0207
##    720       50.4559         67.6265     0.0100   -0.0792
##    740       49.9923         68.1120     0.0100    0.0101
##    760       49.6616         68.2458     0.0100    0.0180
##    780       49.4424         68.4961     0.0100   -0.0425
##    800       49.0898         68.6180     0.0100    0.0125
##    820       48.6806         68.4983     0.0100   -0.0262
##    840       48.4303         68.7273     0.0100   -0.0407
##    860       48.0850         69.3697     0.0100   -0.0104
##    880       47.6589         69.4833     0.0100   -0.0119
##    900       47.3868         69.9136     0.0100   -0.0164
##    920       47.0484         70.5000     0.0100   -0.0194
##    940       46.6978         70.2430     0.0100   -0.0012
##    960       46.5355         70.8139     0.0100   -0.0087
##    980       46.1201         70.6895     0.0100   -0.0211
##   1000       45.6229         70.8181     0.0100   -0.0228
pred_gbmtrain=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[train,]),
                        obs=target_train)
## Using 166 trees...
pred_gbmtest=data.frame(pred=predict(gbm_model,
                                     newdata = data_model[-train,]),
                        obs=target_test)
## Using 166 trees...
print(paste("gbm test:",Gini(pred_gbmtest)))
## [1] "gbm test: 0.845622319906236"
print(paste("gbm train:",Gini(pred_gbmtrain)))
## [1] "gbm train: 0.861493418882915"
print(paste("glm:",Gini(pred_glmnettest)))
## [1] "glm: 0.591052172672043"
plot(gbm_model,i.var="TP60AGE114")

Pour en savoir plus sur la simplification des modèles de machine learning complexes, voici un article sur “interpretable machine learning” IML

Dans le cas spécifique des GBM des méthodes adhoc existent parce qu’un arbre peut être vu comme un ensemble d’indicatrices sur des croisements de variables.
Il y a potentiellement autant de croisements que la profondeur de l’arbre.
Si on décide d’utiliser des “stumps” arbres de profondeur 1 alors chaque arbre est une indicatrice sur 1 variable. Autrement dit ce sont des GLM (parce qu’il y a quand même une fonction de lien) avec des interactions d’ordre 1.
Si on décide d’utiliser des arbres de profondeur 2 alors les arbres construisent des indicatrices sur des croisements de 2 variables. Autrement dit ce sont des GLM avec des interactions d’ordre 2.
Or un GBM (idem pr random forest) est une somme pondérée d’arbres de décisions.
Et une somme pondérée (combinaison linéaire) de GLM est encore un GLM.
Bref on peut ramener un GBM de profondeur 1 ou 2 à un GLM avec des interactions d’ordre 1 ou 2.
Ceci dit il y aura beaucoup de coefficients donc on est plus proche d’une GAM non paramétrique que d’un GLM (au lieu de splines on utilise des fonctions constantes par morceaux)